HIGHLY ANISOTROPIC TEMPERATURE BALANCE EQUATION AND 
ITS ASYMPTOTIC-PRESERVING RESOLUTION. 

ALEXEI LOZINSKI, JACEK NARSKI, CLAUDIA NEGULESCU 

Abstract. This paper deals with the numerical study of a nonlinear, strongly anisotropic 
CN heat equation. The use of standard schemes in this situation leads to poor results, due 

to the high anisotropy. An Asymptotic-Preserving method is introduced in this paper, 
which is second-order accurate in both, temporal and spacial variables. The discretization 
in time is done using an L-stable Runge-Kutta scheme. The convergence of the method 
is shown to be independent of the anisotropy parameter < e < 1, and this for fixed 
coarse Cartesian grids and for variable anisotropy directions. The context of this work 
are magnetically confined fusion plasmas. 
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1. Introduction 

Magnetically confined plasmas are characterized by highly anisotropic properties induced 
by the applied strong magnetic field. Indeed, the charged particles constituting the plasma 
move rapidly around the magnetic field lines, their transverse motion away from the field 
lines being constrained by the Lorentz force. In contrast, their motion along the field 
lines is relatively unconstrained, so that rather rapid dynamics along the magnetic fields 
occurs. This results in an extremely large ratio of the parallel to the transverse thermal 



conductivities, as well as of other parameters characterizing the plasma evolution. 



A prototype simplified model for the heat diffusion in a magnetically confined plasma 



-^ can be expressed by the following nonlinear, degenerate parabolic equation 

^ ^w-V r («n(«)V||w)-V ± -(« ± V ± u) = 0, (1.1) 

where the subscripts || (resp. _L) refer to the direction parallel (resp. perpendicular) to the 
magnetic field lines and u designates the temperature. In writing out the equation above we 
have ignored some important physical phenomena coming from convection and turbulence. 
Nevertheless, our equation contains some important features inherited from the full model 
that lead to substantial difficulties in the numerical treatment of both the full model and 
our simplified one. The diffusion in the direction perpendicular to the magnetic field lines is 
usually dominated by the anomalous transport and the corresponding coefficient K± can be 
taken temperature independent. On the other hand, the coefficient describing the diffusion 
in the direction parallel to the magnetic field lines, k\\, is normally much larger and strongly 
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temperature dependent. It can be described by the Spitzer-Harm law K\\(u) = kqvPI 2 [2^] . 
Moreover, plasma temperatures are extremely high, so that this diffusion coefficient can 
become very big. Passing to non-dimensional variables, we shall write therefore the law 
for K\\ as 

«||(«) = -M 5/2 , 

where e is a small parameter, < £ < 1. An accurate resolution of the parallel and 
perpendicular diffusion processes plays a crucial role in understanding of the plasma dy- 
namics and the energy transport phenomena. It is therefore very important to develop 



and to study efficient numerical schemes to solve problem (1.1). It is also desirable to 
have a scheme that works robustly for all values of e from e << 1 to e ^ 0(1) since this 
parameter enters the equation in combination with a non-linear term so that the effective 
value of the diffusion coefficient can vary strongly over the computation domain following 
the variations in u. This is the primary motivation of the present work. 



Anisotropic, nonlinear diffusion equations of the type (1.1) arise in several other fields of 
application and a lot of efforts were made to construct efficient numerical methods for this 
challenging problem. To mention some examples, such non-linear evolution equations of 
parabolic type occur in the description of isentropic gas flows through a porous media [2] 
or in the description of transport phenomena in heterogeneous geologic formations, such as 
fractured rock systems [5] , which are of fundamental interest for petroleum or groundwater 
engineering. In addition, these equations appear also in image processing, related to the 
elimination of noise and small-scale details from an image [21 Q21 |23] or in the description 
of the anisotropic water diffusion in tissues of the nervous system [I] . 



From a numerical point of view, problems of the type (1.1) are very challenging, as one 
deals with singularly perturbed problems, the model changing its type in the limit e — > 0. 
Standard schemes suffer from the presence of very ill conditionned matrices (typically with 
a condition number of order I /(eh 2 ) where h is the discretization step in space). Solving 
an equation with such a matrix on a computer accumulates the rounding errors and may 
lead to completely wrong results. Note that this drawback cannot be overcome by a mesh 
refinement since it results only in worsening the condition numbers of the matrices in the 
discretized problem. 

Several methods were investigated in literature to cope with this type of anisotropic 
problems, using for example high order finite element schemes [TO], preconditioned conju- 
gate gradient methods in a mixed spectral/finite difference scheme [15] or introducing an 
artificial "sound" method, to represent the fast thermal equilibrium along the field lines 
|17j . All these methods however are rather involved and moreover their range of appli- 
cation is limited, as they are efficient only until a threshold value for e, and cannot thus 
recover the limit regime e — > 0. Another class of employed numerical methods are hybrid 
strategies, which consist in coupling different numerical schemes valid in different regions 
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of the domain. For example in this case, one can couple the resolution of the singular 
perturbation problem there where e ~ 0(1) with the resolution of a limit problem for 
e -^ 1. These methods suffer however from the fact, that the coupling conditions between 
the two models are hard to establish and the interface between the two regions difficult to 
localize. 

The objective of the present paper is to introduce an efficient numerical scheme based 
on the Asymptotic-Preserving methodology, which allows for an accurate resolution of the 
singularly perturbed problem, uniformly in e, with little additional computational cost, 
and using a grid which is not necessarily aligned with the magnetic field, so that one can 
exploit simple Cartesian grids, for example. Initially, AP-techniques were introduced in 
[12] . to deal with singularly perturbed kinetic models. The key idea is to reformulate the 
singularly perturbed problem into an equivalent problem, which is however well-posed if 
we set e = there. The reformulation of the here proposed method is based, similarly as in 
rrp] , on introducing a new auxiliary variable, as proposed earlier in an elliptic framework 
in [7j, and replacing the terms of the equation multiplied by 1/e are by the new terms 
with an 0(1) factor. From a numerical point of view, this procedure means transforming 
the problem of condition number ~ 1/e, into a well-conditioned problem, which switches 
automatically from the singularly perturbed problem to the limit problem, as e — > 0. 

The difference between the method presented in [16] lies first in the treatment of the 
non-linearity. Instead of fixed point iterations used to approach the nonlinearity (w n+1 ) 5//2 , 



we choose to implement a much simpler linear extrapolation method (see Section |3.2 for 



details). Moreover, we develop here a robust asymptotic-preserving scheme of second order 
in time, which has no analogue in the existing literature, to the best of our knowledge. 
Finally, this paper contains also a detailed mathematical study of the problem. 

The paper is organized as follows: Section [2] contains a description of the problem com- 
pleted by a mathematical study. In Section |3j we present the numerical method based on an 
asymptotic preserving space discretization and develop three different time-discretizations: 
implicit Euler, Crank-Nicolson and L-stable Runge-Kutta methods. Finally, in Section [4] 
we present some numerical results, focusing on the AP-property of the schemes. 



2. Description of the problem and mathematical study 

We consider a two or three dimensional anisotropic, nonlinear heat problem, given on 
a sufficiently smooth, bounded domain Q C M d , d = 2, 3 with boundary I\ The direction 
of the anisotropy is defined by the time-independent vector field b G (C OQ (Q)) d , satisfying 
\b(x)\ = 1 for all igfl. 

Given this vector field b, one can decompose now vectors v G M d , gradients V0, with (j)(x) 
a scalar function, and divergences V • v, with v(x) a vector field, into a part parallel to the 
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anisotropy direction and a part perpendicular to it. These parts are defined as follows: 

v || := (v ■ b)b , v± := (Id — b (g> 6)t> , such that f = v\\ + v j_ , 

V||0 := (b- V0)6, V ± := (Id - b ®b)V(p , such that V0 = V\\(j) + V ± 0, 

V|| • v := V • f || , Vj_ • v :— V • v± , such that V • v = V\\ ■ v + Vj_ • v , 

(2.2) 
where we denoted by (g> the vector tensor product. 

The boundary T can be decomposed into three components following the sign of the 
intersection with b: 

T|| := {x e T / b(x) -n(x) = 0}, 

v in ■= {x e r / b(x) ■ n(x) < 0} , r out ■.= {x e r / 6(x) • n(x) > 0} , 

and I\ = T in U T out . The vector n is here the unit outward normal on F. 

With these notations we can now introduce the mathematical problem, we are interested 
to study. We are searching for the particle (ions or electrons) temperature u(t, x), solution 
of the evolution equation 

' d t u - JV|| ■ (A||U 5 / 2 V||u) - V± ■ (A± V±u) = , in [0,7] x to, 

\n\\ ■ 0V 5 / 2 (£, -)V||u(*, •)) + n± ■ (A ± V±u(t, •)) = ~iu{t, ■) , on [0, 7} x V ± , 

V ± u(t,-) = 0, on [0,T]xr,|, 

«(0, •) = «°(0 , in fl. 

(2.3) 
The coefficient 7 is zero for electrons and 7 > for ions [211 121] • The problem (2.3) 
describes the diffusion of an initial temperature u° within the time interval [0, 7] and its 
outflow through the boundary I\. Let us denote in the following the time-space cylinder 
by Qt '■= (0, 7) x Q. The parameter < e ^C 1 can be very small and is responsible for the 
high anisotropy of the problem. We shall suppose all along this paper, that the coefficients 
A\\ and A± are of the same order of magnitude, satisfying 

Hypothesis 1. Let T± consist of two connected components Ti n = {x G T/n ■ b < 0} and 

Font = {x E T/n ■ b > 0} such that 

either 

case A: n = —b on Y in and n = b on T out 
or 

case B: n ■ b > —k on T in and n ■ b < k on T out with some constant < x < 1. 
All the components T in , T out and Fu are sufficiently smooth. We suppose moreover < s < 
1 and 7 > fixed. 7he diffusion coefficients An G W 1 ' 00 ^) and A± G M dxd (l^ 1 ' 00 (n)) are 



(P) 
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supposed to satisfy 



< A < A\i(x) < A 1 , f.a.a. ieO, 
A)|M| 2 < v t A±(x)v < v4i||t;|| 2 , Vw G M d and f.a.a. x e Q, 



(2.4) 
(2.5) 



with < A < A\ some constants. 



Putting formally e = in (2.3) leads to the following ill-posed problem, admitting 
infinitely many solutions 



' -V r (A||M 5 / 2 V||M) = 0, in [0,T]xQ, 
n r (A ll u 5 / 2 (t,-)V ll u(t,-))=0, on [0, T] x T ± 
V_lu(V) = 0, on [0,T]xr,|, 

k u(0,-) =«°(-), in tt. 



(2.6) 



Indeed, all functions which are constant along the field lines, meaning VmU = 0, and 
satisfying moreover the boundary condition on r\[, are solutions of this problem. From a 
numerical point of view, this ill-posedness in the limit e — > can be detected by the fact, 



that trying to solve (2.3) with standard schemes leads to a linear system, which is very 



ill-conditioned for < e <C 1, in particular with a condition number of the order of 1/e. 
The aim of this paper will be to introduce an efficient numerical method, permitting 



to solve (2.3) accurately on a coarse Cartesian grid, which has not to be adapted to the 
field lines of b and whose mesh size is independent of the value of e. The here proposed 
scheme belongs to the category of Asymptotic-Preserving schemes, meaning they are stable 
independently of the small parameter e and consistent with the limit problem, if e tends 
to zero. The construction of the here developed AP-scheme is an adaptation of a method 
introduced by the authors in an elliptic framework (see [7j), to the here considered non- 
linear and time-dependent problem, and is based on a reformulation of the singularly 



perturbed problem (2.3) into an equivalent problem, which appears to be well-posed in 
the limit e — > 0. But before introducing the AP-approach, we will start by studying in 



the following section the mathematical properties of problem (2.3). The test configuration 



chosen all along this paper is the diffusion of an initial temperature hot spot (see Figure 
fl| along arbitrary magnetic field lines b. 



2.1. Mathematical properties. Before starting with the presentation of the AP numer- 



ical scheme, let us first check some properties of the diffusion problem (2.3) for fixed e > 0. 
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Figure 1. Diffusion of a hot temperature spot along the magnetic field lines. 



For notational simplicity, we consider a slightly more general form of problem (P) 



(Prr 



' d t u - V|| ■ (A\\\u\ m - l V\\u) - Vl • (A ± V ± u) = , in [0,7] x to, 

Au\u\ m - 1 ni r V\\u + A ± n ± -V±u = ->yu J on [0,T]xr_L, 
V±u = 0, on [0,T]xr N , 



(2.7) 



{ «(0, 



u 



in fi . 



by setting m = 5/2+1 and redefining A\\ 
as -An for any e > 0. Equations of the type (2.7) are rather well studied in the literature. 



for any m > 1. We obtain the particular case (2.3 



We refer to the classical works [HJ [9j [13] as well as to the more modern literature on "The 
porous medium equation" as reviewed in (TJ [22] . However, all these references normally 
treat only an isotropic version of the problem above, i.e. the non-linearity of the type 
is present in front of all the derivatives of u. An anisotropic equation of the form 



u 



m—l 



(2.7) is studied in [TT], but only in the case m < |^-, so that the value of m pertinent 
to our application is not covered. Another feature of our setting, which is not sufficiently 
covered in the existing literature, is the prescription of Robin boundary conditions. This 
is the reason why we wish to study in this paper the existence, uniqueness and positivity 



of solutions to (2.7). 



We shall first introduce the concept of weak solution of problem (2.7) and state the 



existence/uniqueness theorem. Note that unlike the literature cited above, we assume 
from the beginning that the initial conditions are bounded and strictly positive, and prove 
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the same properties for the weak solutions. Our treatment is thus performed under much 
less general assumptions than usually required, but this is quite enough for our application. 

Definition 2.1. (Weak solution) Let u° e L°°(fi). Then ueW with 

W :={uE £°°(<2oo), such that VT > 

V ± u G L 2 (Q T ) , \u\ m ~ l V\\u G L 2 (Q T ) , 8 t u G L 2 (0, T; (H\Q))*)} , 



is called a weak solution of problem (2.7), ifu(0, •) = u° and if for all T > one has 

/ (d t u(t,-),(f)(t,-)) iH i)* t mdt+ / A|||M| m_1 V||M- V\\(f)dxdt (2.1 

Jo Jo Jn 

+ / A ± V±u-V ± (j)dxdt + >y / ucj)dadt = 0, V0 G X> 



Remark 2.2. AZ/ the terms in this variational formulation are well-defined for any u EW 
and cj) G V. Indeed, as u G L°°(Qt), Vj_m G L 2 (Q t ) and |M| m_1 V||M G L 2 (Q T ), one also 
has \u\ m ~ l Vu G L 2 (Q T ), and thus \u\ m G L 2 (0,T; H\n)). This means that \u\ m has a 
trace on dVl, belonging to L 2 (dVt). As L 2m (dVt) C L 2 (dVt) for all m > 1 (D, is bounded), 



one has u\gn G L 2 (dQ), justifying thus the boundary integral in (2.8). Moreover, we have 



the continuous inclusion W C C([0, T]; (H 1 ^))*). This shows that one can impose the 
initial condition w(0, •) = u° with u° G L°°(fi) C (H 1 ^))*. 

Remark 2.3. Actually, we have a sharper characterization of continuity in time for func- 
tions in W. Indeed, let Y be a Banach space and let us denote by C([0, T]; Y w ) the space of 
weakly continuous functions with values in Y , which means for each ip G Y* the mapping 
t i — > {ip,u(t))Y*,Y is continuous. Then, one can prove (see [2] for more details) that 

WcL-(0,T;L 2 (fi))nC([0,T];(if 1 (fi)r)cC([0,r];L 2 y (fi)). 

Theorem 2.4. (Existence / 'Uniqueness /Positivity) Let u° G L°°(fi) satisfy < (3 < 

u° < M < oo on Vt, for some (3 > 0. Assume Hypothesis [7] and m > 1. Then there exists 



a unique weak solution u G W of (2.1), which satisfies ce < u < M a.e. on Q^, with 



some sufficiently small c > and some sufficiently large K > 0. 

Before proving this theorem, let us also define the sub- and super-solutions to problem 



(2.7) and establish a comparison principle for them. 



Definition 2.5. (Sub/ super- solutions) A function u G W is called a weak sub- (resp. 



super-) solution to problem (2.1) if the variational formulation (2.8) is verified for allcj) ET> 



with (f) > on Qoo, and where the equality sign is replaced by < (resp. >). 
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Lemma 1. (Comparison principle) Assume Hypothesis [7] and m > 1. Lei «i be 



a non-negative sub-solution and u 2 be a non-negative super- solution to (2.1) such that 
i*i (0, x) < 1*2(0, x) /or a. a. a; G fi. // at least one of the functions U\, 1*2 is strictly positive, 
i.e. VT > 3pV > such that U\ > pV or u 2 > 0t on Qt, then U\ < u 2 on Q co- 
Proof. For any k > 0, introduce the function H k : R — > M. as 

f 0, if u < 
Hkiu) = I ku, if < u < I 
I 1, ifu>£ 

and put = F fe (wf - it^). Note that G L 2 (0, T; H l {Vt)) since i*i,ii 2 G W and thus 
Vw™, Vm™ G L 2 (Q T ). Observe also that the gradient of <fi is zero outside from the set 
u t = {(ti x ) e Qt '■ < u™ — u™ < r}- Choosing this as the test function in the 



inequalities (2.8) for u\ and 1x2 and subtracting the second one form the first one gives 
/ (dt(ui - u 2 ), H k (u™ - v%)) {m y, H i dt < -A [J A|||V||« - u™)\ 2 dxdt (2.9) 



-mk A ± V x (u x - u 2 ) ■ [«i Vi(«i - 1*2) + {u™- 1 - u r 2 n - l )V ± u 2 }dxdt 



,k 



-7/ f ( Ul -u 2 )H k (u^-u^)dadt 
Jo ir. 



<mk A 1 _\V±{u 1 - 1*2) • V ± w 2 |« - u™~ l )dxdt 






T 



since (u± — u 2 ) and Hk(u™ — u™) are of the same sign. We have moreover 

„.m _ „.m p 

v-i u 2 ^ o m ■ ^ 

Mx + M 2 ^P 

on u)\ with a constant C m > depending only on m. Indeed, the first inequlity here holds 
for any U\ > u 2 > and the second inequality follows by noting that u™ — u™ < l/k on 



u^ and U\ +u 2 > j3 T on Q T . We see that (2.9) combined with the inequality above implies 
J (d t {u 1 -u 2 ),H k (uT-u 2 n )) {H iy >H idt< r ^ II A^V^-^-V^ldxdt. 



Let us now take the limit k — > 00 in this inequality. We have meas(u}^) — > so that 

limsup / (a i (n 1 -ii 2 ),^«-«™))(Hi)*,Hirft<0. (2.10) 

fc— s>oo Jo 
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On the other hand, 



lim 



k~ >oo 



(d t ( Ul -u 2 ),H k (u?-u™)) {H 



(Hi)*,H 



dt = (ui—u 2 ) + (T,x)dx— I {u\— u 2 ) + (0,x)dx. 
Jn Jn 

(2.11) 

where (u) + = (u+\u\)/2 denotes the positive part of u. Indeed, Hk{u™ — u™) — > H{u\ — u 2 ) 

a.a. on Qt where H denotes the Heaviside function {H(x) = 1 for x > and H (x) = for 

x < 0). Observing that d t {u\ — u 2 ) + = d t {ui — u 2 )H{u\ — u 2 ) in the sense of distributions, 



proves (2.11) for sufficiently smooth U\,u 2 . A standard density argument shows then 



that (2.11) actually holds for any Ui,u 2 G W. Note, in particular, that the terms at 



the right-hand side of (2.11) are well defined for functions in W thanks to the inclusion 



W C C([0,T];Ll {ft)), cf. Remark [2T3J Comparing ( ]2.10[ ) and fl2.11| ) and taking into 
account {ui — u 2 ) + = on Q at t = 0, yields 

{u 1 -u 2 ) + (T,x)dx<0, (2.12) 

which implies U\ < u 2 on Qoo- D 

The construction of the following remarkable sub-solution is essentially due to M. Pierre 



Lemma 2. (Construction of a weak solution) Assume Hypothesis [7] and m > 1. For 

any (3 > 0, there exists a weak sub-solution w to problem (2.7) satisfying c < w(0,x) < (3 
for x G Q and w(t, x) > ce~ Kt for (t, x) G Qoo, with some constants c,K>0 which depend 
only on (3. 

Proof. We will construct a smooth sub-solution w satisfying all the announced properties. 
We thus rewrite the definition of a sub-solution in the strong form supposing from the 
beginning that w > 0: 



d t w Vii ■ {AnVnw™) - V ± • {A ± V ± w) < 0, 

m 

— A\\n\\ ■ V\\w m + A ± n ± ■ V±w + jw < 0, 

fit 

n± ■ Vj_u> < 0, 



on 



(0, oo) x n 

(0,oo)xr ± , 
on (0, oo) x T|| 



on 



(2.13) 

(2.14) 
(2.15) 



The construction of such a function w will be performed separately for the two cases 
mentioned in Hypothesis [TJ 

Case A: One can construct in this case a new coordinate system £i, . . . ,£d on Q such 
that the coordinate lines ^ coincide with the 6-field lines and the surfaces £$ = const are 
perpendicular to these lines. Domain Q is represented in these coordinates by a cylinder 
Cj ? = D x (0, 1) with f = (&, . . . , ^_i) G D and £ d G (0, 1). We thus have V|| = bx^ with 
some scalar strictly positive field x- We assume that the component T in of the boundary 
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is represented by D x {£ d = 0}, T out is represented by D x {£ d = 1} and T\\ is represented 
by 3D x (0,1). 

We are searching now for a sub-solution under the form w(t, x) = <5(£)(sin(7r£ rf ) +r](t)) 1 ^ m 
where 5(t) and r)(t) are two positive decreasing functions which are yet to be adjusted. We 
observe immediately that V±w = on Q for all time so that (2.15) is automatically 
satisfied. The remaining boundary conditions (2.14) should be checked on Ti n and T out . 
We remind that n = n\\ = b on Y out (t; d 
Substituting the Ansatz for w into (2.14) now gives 

1 



1). Similarly, n 



mi 



0). 



m 



-An X 5 m n + 7077™ < 0, for £ d = and £ d = 1. 



This holds if one takes i] = Ki8 m ^ n ~' 1 ' where K\ = (minggjv ^-A\\ X ) > 0. 

It remains to check (2.13). Substituting the Ansatz for w, this inequality is reduced to 

d 



5(sm(ir£ d ) + r])< 



5 m d 5 m 

— X^" (-A[[X) 7rcos« d ) H A||X 2 7r 2 sin(7r£ d ) < 0. 

m ot, d v 7 m 



(2.16) 



Note that we have denoted the time derivative here by a dot and we neglected a term 
with r) since it is negative (the function i](t) is decreasing). We divide now both sides by 
(sin(7r£d) + r})™ and bound each term on the left-hand side as 



S m 7rcos(^ d )x d(A l]X ) < tt5^ X 



m(sin(7T^) +r])m d£, 
and 



mr/' 



d(A llX ) 



% 



ir5 m x 



mK^S 



m fm-1 



d(A ]lX ) 



d£a 



5- 



IX 



(? 

m 



sin(7r£ d 



< 



(V 



sin« d ) +T])< 



in 



A llX 2 7T 2 (sm(n^ d )) 



i-- 



(V 



mKl 



2„2 



— X 



d(A\\ X ) 



d£ a 



< —A llX z n 
m 



We see now that inequality (2.16) will be satisfied if we require 

S + K 2 5 + K 3 5 m < 
with 



(2.17) 



Ko 



71 



max 



x k {A,,x) 



and K* 



n 



One can thus take 5(t) = 5 e" {K2+Ki)t with any 5 G (0, 1]. 



max L4ii v 
m £en e ' M 



m(m-l) - m (m-\)(Ki+K$)t\ 



is a sub- 



In summary, w(t,x) = 5$e ^ 2+jft:3 ' )t (sin(7r^) + K\a t 
solution. Clearly, for any [3 > one can take 5 small enough so that w(0,x) < (3. 



Moreover, for any t, w(t,x) > S n K 1 < m e - m ( K 2+K 3 )t go ^^ we j^yg p r0 ved the statement 
of the Lemma putting c = S n K 1 , if = m(K 2 + K 5 ). 

Case B: Let G C 2 (f2) be a strictly positive function such that <f)(x) > 1 on 0, n± ■ Vj_0 = 
n ■ V0 < on T|| and V0 = — C<M on T_|_ with some sufficiently big constant ( > 0, to be 
prescribed later. 
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We are searching now for a sub-solution under the form w(t, x) = ce~ Kt <j)(x) where c, K 
are some positive constants which are yet to be adjusted. We observe immediately that 
(2.15) is automatically satisfied for such w. The left-hand side of (2.14) can be written as 

rm -e- mKt A\\n\\ ■ V,|0 m + ce- Kt (A ± n ± ■ V ± + 7 0) 



m 



< ce~ Kt (-Cc m -V (m - 1)Kt ,4|| \n\\ \ 2 <p m - (A ± n ± ■ n ± + 7 
<ce- Kt (-(A ± n ± -n ± + -f<j)) 

and thus it is negative provided ( is chosen sufficiently big. Indeed, A±n± ■ n± is uniformly 
bounded from below by a positive constant in view of the geometrical hypothesis of case 
B. 



It remains to check (2.13). Substituting the Ansatz for w into this inequality yields 

/ m-l -(m-l)Kt \ 

ce- Kt l-Kcf> V|| ■ (A||V||<r) - Vi. ■ (A ± V ± 0) J < 0. (2.18) 

This inequality is satisfied provided we take c < 1 and 

K = — maxlVii • (AnVn^ m ))\ + max|V± • (A ± V±(p)\. 

Finally, for any /3 > one can take c small enough so that w(0, x) = c<f)(x) < (3. Lemma 
is thus proved also in case B. 

□ 

Remark 2.6. In the case of a simple "aligned" geometry, i.e. b = Cd and Q = Dx]0,L[ 
with D a domain in IR d-1 ; and supposing An = const, one can easily construct a sub- 
solution satisfying a sharper estimate: under the assumptions of the preceding Lemma, 
there is a sub-solution such that 

w[t,x) > s— ■ 

K J ~ (1 + Kt)^ 

Indeed, one can repeat the proof as in case A of the preceding Lemma, taking £' = (xi, . . . , Xd-i), 
£,d = Xd/L, up to the differential inequality (2.11). One observes now that K2 = so that 
one can take 8(t) = — 1 with any 5q > and K = (m — l)K 3 5^ 1 . Our sub-solution 

j_ 
A ( Ki5 ml ' m ~ 1 ' l \ m K™& m 

is thus w = ° 1 sm(iiXd) H — n ?^\m and w > l — l ±m— as stated. 

(l + Kt)m=T ^ {l+Kt) J (l+Xt)^=T 

Let us now turn to the proof of our main result. 



Proof of Theorem 2-4- We shall first regularize the problem, in order to avoid the 
degeneracy. Then, in a second step, we shall treat the non-linearity via a fixed point 
argument. Finally, a priori estimates shall help us to pass to the limit in the regularized 
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problem, to prove existence. The comparison principle above will be used to establish the 
uniqueness and the positivity of the solution. Let us now detail these steps. 

1st step: Regularization 
Fix < a < 1 and assume that M > is an upper bound for u°. Introduce for notational 
simplicity the following functions a a ,A a : R — > M. + 

a a (u) := [a + min(|w|, M)] m_1 , A a (u) :— / a a (s) ds , 

Jo 

and consider the regularized version of (2.8): find u a G W\{$$,T\ H X (V£), L 2 (fl)) (i.e. u a G 
L 2 (0,T;H l (Q)) and d t u a G L 2 (0,T; (H 1 ^))*)) such that u a {0, -) = u° and 

(d t u a {t,-),(f)(t,-)) iH i ) * tH idt+ / / A\\a a (u a )V \\u a ■ V\\<pdxdt 
T Jo Jn t (219) 

+ / A ± \7 ±u a ■ V ±(j) dxdt + ^ / u a (frdadt = 0, V0 G £>. 

Jo Jn Jo Jr x 

By standard arguments, this problem is well posed. Indeed, consider the mapping 

T:B R (0) -»- 5 fl (0), S fl (0) := {^ G L 2 (Q T ) / \\v\\ l2{Qt) < R} , 

where we associate to any v G -Br(O) the unique solution u G ^^(OjT; H 1 ^), L 2 (Q)) of 
the linearized, regular parabolic problem 

/ (d t u(t,-),(f)(t,-))( H i)* iH idt+ A\\a a {v)V\\u-V\\(f)dxdt 

Jo Jo Jn 

+ A ± V ±u ■ V ±(f) dxdt + 7/ / u(t>dadt = 0, \/(peV. 

Jo Jn Jo Jr ± 



Indeed, taking R := a/T||m || 2 , the mapping T is well-defined, continuous and T(B R (0)) 
is relatively compact in L 2 (Qt)- The continuity follows from the fact that for v n — > v 
in L 2 (Q T ) and v n ^ w in W^iO, T;i7 1 (r2),L 2 (r2)), the Lebesgue dominated convergence 
theorem permits us to pass to the limit in the linearized term of the variational formulation. 
By Schauder fixed point theorem, T has a fixed point T(u) = u, which provides a solution 
to pi*!. 



The solution u a of problem (2.19) satisfies < u a < M, provided we have < u° < M. 
Indeed, define u~ := min(0,M a ) < 0. Then one gets for the initial condition u~(0,-) = 0. 
Taking u~ as the test function in the variational formulation (2.19), yields immediately 

/ Ur(T.x)\ 2 d,x + / 



- / \u a (T, x)\ 2 dx + A\\a a (u a )\V\\u a \ 2 dxdt 

2 Jn Jo Jn 

+ / A ± \V±u~\ 2 dxdt + 'j / \u~\ 2 dadr = 

Jo Jn Jo Jr± 
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which shows that u~(T, ■) = 0. Since the same argument can be applied to any final time 
T, we have u a > in Qoo- 

To prove the estimate from above, define w+ := max(0, u a — M). Observe that u+ (0, ■) = 
and take u^ as the test function in the variational formulation (2.19): 

1 



|w„(T, x)\ dx + 



A\\a a (ua)\V\\Ua\ dxdt 



o Jn 

T 



- I / Aj_|V ±u+r dxdt + 7 

'o Jn Jo Jr A 



u a u~a da dr = 0, 



which shows that u^(T, ■) = 0. Since again the same argument can be applied to any final 
time, we have u a < M in Q^. 
2nd step: A priori estimates 
In order to pass to the limit a — > 0, we will need some a priori estimates for the solution 

i the variational formu 

A\\a n (u n )W \\u n \ 2 dxdt 



u a , independent of a. Taking <p = u a in the variational formulation (2.19) yields 
\u a (T, x)\ 2 dx + 



" I \u a [;i',x)\*dxi- I I Alio,,,..,,, 
* Jn Jo Jn 

+ / A±\V ±u a \ 2 dxdt + 7 
Jo Jn Jo Jr 



\u a \ 2 da dt = - / \u°(x)\ 2 dx, 
2 Jn 



(2.20) 



which implies 

IMIl°°(o,t ; l 2 (q)) < ||m°||l 2 (q), J J Q a a (u a )\V\\u a \ 2 dxdt < C||w°||^2 (n) 

||V_lM q ||l2(q t ) < C||w°||x^(fj), ||Wa|U 2 ([0,T]xr x ) < C|| u °l|L 2 (n), 

with some constant C > 0. 

Taking now = A a (u a ) in ( |2.19| ), which is permitted since w a e L°°(<2 T )nL 2 (0,T; H 1 ^)), 
yields 

(dtU a ,A a (u a )}(Hi)*,mdt+ \ \ An\V\\(A a (u a ))\ 2 dxdt 



o jo Jn 

+ / / A±a a (u a )\V±u a \ 2 dxdt + 7 / / u a A a (u a )da dr = . 

The first term can be rewritten as 

/ (d t u a ,Aa(u a ))(Hi)* tH idt= / $? a (u a (T,x))dx- / ty a (u°(x))dx , 

Jo Jn Jn 

with ty a {u) '■= Jq A a (s) ds. Due to the facts that < u a < M, A a (u a ) > and ^/ a (u a ) > 0, 
we get 

/ / A\\\V\\{A a {u a ))\ 2 dxdt < C [ V a (u°(x))dx< ° I \a + u°(x)) m+1 dx 

Jo Jn Jn m{m + 1) J n 

< C (a m+1 + M m+l ) . 
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Thus, we have that the family {V\\A a (u a )} a , a g]0, 1] is bounded in L 2 {Q T ). More- 
over, {V±A a (u a )} a is also bounded in L 2 (Q T ), since V±(A a (u a )) = a a (u a )'V±u a and 
a a{ u a) is uniformly bounded by some positive constant. Hence, {A a (u a )} a is bounded in 
L 2 (0,T;H\Q)). 

Let V = H 1 (fl)nL°°(fi) be the Banach space with the norm ||-||y = ||- ||h 1 (q) + ||'||l°°(q)- 
For any in L°°(0,T;V), 



T 



{d t A a (u a ),4>)( H i)* <H idt 



(d t u a , a a {u a )4>){mY,m dt 



< C\\(/)\\l°°(0,T;V) 



with a constant C independent of a. This follows from the variational formulation (2.19) 
with replaced by a a (u a )(j) and from the estimates (2.20). We see thus that the family 
{d t A a (u a )}a is bounded in L 1 (0,T;l / *). We remind also that {A a (u a )} a is bounded in 
L 2 (0,T; if 1 (fi)). Aubin-Simon compactness lemma [20] applied to the triple of spaces 
H 1 ^) C L 2 (Q) C V* allows us now to conclude that the set {A a (u a )} a is relatively 
compact in L 2 (0,T; L 2 (Q)) = L 2 (Q T ). 

3rd step: Passage to the limit 
The aim now is to pass to the limit a — > in the variational formulation (2.19) in order to 
show the existence of a weak solution of problem (2.3). The a priori estimates of the last 
step permit us to show, that there exists a function u G L 2 (Q T ) satisfying < u < M in 
Qt and such that after extracting a sub-sequence from {u a } a , we have 

u a ^a^o u in L 2 (Q T ), u a \ r± ^ a -+ u\ r± in L 2 ([0,T] x T±) , 

V_l« q ^ a ^o V_lu in L 2 (Q T ) and d t u a — V->o d t u in L 2 (0,T; (H 1 ^))*) . 

To pass to the limit in the non-linear term, we use first the fact that {A a (u a )} a is bounded 
in L 2 (0,T; if 1 (r2)), so that there exists some function w G L 2 (0,T; H l (Q)) satisfying 

A n (« a )^o» m L 2 (0,T;H\Q)). 

In order to identify the function w, we need some pointwise convergence of the sequence 
u a - For this, recall that the sequence {A a (u a )} a is relatively compact in L 2 {Qt). Thus, 
up to a sub-sequence 

A a (u a ) — j-q.^o w i n L 2 (Q T ) , hence A a {u a ) — > a ^o w , a - e - i n Qt • 



-^ Q ^o A l w 



u 



Since for any fixed u G [0,M], A a («) — t- q ^o A{u) := — M m , we have u a 

a.e. in Qt- This permits to identify the function A -1 w; with u, so that w — — 

L 2 (0,T;H\n)). 

All these convergences allow now to pass to the limit in the variational formulation (2.19) 
and to conclude the proof of existence for problem (2.7). 

4th step: Positivity and uniqueness 
Let u be a weak solution to (2.7). We use first the sub-solution constructed in Lemma 
[2] and the comparison principle in Lemma [I] to verify that u > ce~ Kt with some positive 
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constants c small enough and K large enough. We remark then that M is a super-solution 
to (2.7). Again, by Lemma IT] we see that u < M. 

Suppose now that ( 2.8[ ) admits two solutions U\ and u 2 in W with the same initial 
condition U\ = u 2 = u° at t = 0. We know already that they are both strictly positive, 
so that Lemma [T] implies U\ < u 2 on Q,^. Since the two solutions U\ and u 2 are perfectly 
interchangeable in the above argument, we have also U\ > u 2 and thus U\ = u 2 on Qoo- D 

Let us remark here, that due to the strict positivity of the solution, in particular to the 
property that u > ce~ Kt a.e. on Qoo, with some c > and K > 0, we have 



Corollary 2.7. Under the hypothesis of theorem 2.4 • the weak solution u of (2.1) belongs 
to the more regular space 

W := {u G L°°(<2oo), such that VT > : 116 L 2 (0, T; H\Q)) , ftw G L 2 (0, T; (H 1 ^))*)} . 

Moreover, one has 4;\\u(t, -)\\L 2 {n) < . 

3. Numerical method 

3.1. Semi-discretization in space. The singular perturbation problem ( 2.3[ ) is a highly 
anisotropic equation. Its variational formulation reads: find u(t, ■) G V := H l (Q) such that 

(P) (d t u{t,-) 1 v) v * y + - f A\\\u\ b/2 V\\u{t,-)-V\\vdx (3.21) 



+ A ± V±u(t,-)-V±vdx + >y u(t, -)v da = 0, Vw e V 
in ir ± 

for almost every t E (0,T). As mentioned already in Sectional this problem becomes 
ill-posed if we take formally the limit e — > 0. Indeed, only the leading term survives in this 
limit, so that any function from the space 

Q ■= {p e V / V||p = in Q} 

would be a solution. It is easy to establish, however, the well-posed problem for the limit 
of the solutions to (P) as e — > 0. Indeed, one can restrain the test functions in (P) to be in 
the space Q so that the ^-dependent term disappears and the correct problem in the limit 
e — > reads: find u(t, •) 6 Q such that 

(L) (d t u(t, •), v)v*,v + / A±V±u(t, •) • V±v dx + -f u(t, > da = 0, Vf e Q 

for almost every t G (0, T). 

The discussion above shows that a straight-forward discretization of problem (P) may 
lead to very inaccurate results when e « 1. Indeed, setting e = would result in a singular 
problem, so that the problem with e « 1 would be very ill-conditioned. To cope with this 
difficulty and to obtain a numerical scheme which is uniformly accurate with respect to e, 



16 



A. LOZINSKI, J. NARSKI, C. NEGULESCU 



we introduce an Asymptotic-Preserving reformulation, very similar to the one introduced 



in [7]. The idea is to rewrite the singularly perturbed problem (3.21 ) in an equivalent form, 
which is however well-posed when one sets there formally e = and gives moreover the 
correct limit problem (L). In order to do this, we introduce the auxiliary unknown q by 
the relation eVuq = u 5//2 V||U in f2 and q = on T in , which rescales the nasty part of the 
equation permitting to get rid of the terms of order 0(l/e). The reformulated problem, 
called in the sequel the Asymptotic-Preserving reformulation (AP-problem) reads: find 
(u(t, •), q(t, •)) G V x £, solution of 

(A ± V ±u) ■ V ±v dx + / A||V||g- Vj|Wcfe + 7 



(AP) 



where 



(^,v)v*,v + 



dt 



uvds = 0, 

V^G V 



A\\u b/2 V\\u- V\\wdx-e \ AiiViig • Vnwdx = 0, \/w e C. 



£:={qe L 2 (Q) / Vyg 6 L 2 (Q) and q\ Tin = 0}. 



(3.22) 



(3.23) 



System (3.22) is an equivalent reformulation (for fixed e > 0) of the original P-problem 

in (AP) leads to the well-posed limit problem 



(3.21). Putting now formally e 

,du 
{ ~di 



(L') 



v) 



v*,v 



(A ± \7±u) ■ V±vdx+ / AiiViig- Viiudx + 7 



uvds = 0, 

Vt;G V 



A\\u b/2 V\\u- V\\wdx = 0, VweC, 



(3.24) 

which is equivalent to problem (L). Note that q acts here as a Lagrange multiplier for 
the constraint u G Q, which provides the uniqueness of the solution. Hence the AP- 
reformulation permits a continuous transition from the P-model to the L-model, which 
enables the uniform accuracy of the scheme with respect to e. 

Let us now choose a triangulation of the domain Q with triangles or quadrangles of order 
h and introduce the finite element spaces V/, C V and Ch C C of type P*,. or Q^ on this 



mesh. The finite element discretization of (3.22) writes then: find (uh,qh) G Vh x Ch such 
that 

du h 



Of 



v h dx+ / (A±V±u h ) ■ V±v h dx + 



(AP)h { 



A^WQh -V\\v h dx + 7 
<> Jr 

Vv h G V h 



u h v h ds = 0, 



/ A\\u 5 h /2 \7 )\u h ■ V\\w h dx-s Ai\V\\q h ■ V\\w h dx = 0, \fw G C h . 
Jn Jo. 



(3.25) 
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Remark that this system is continuous in time and also nonlinear, such that one has to 
develop now a procedure for the linearization and the discretization in time. This procedure 
has to be chosen carefully, such that the AP-property developed so far, is not destroyed. 
This is the aim of the next section. 

3.2. Semi-discretization in time. In order to approach numerically the time derivative 



in (3.25), we introduce three different schemes : a standard first order, implicit Euler 
scheme, the Crank-Nicolson scheme and a second order, L-stable Runge-Kutta method. We 
show in the following that the first order Euler-scheme is stable and asymptotic-preserving. 
The Crank-Nicolson schemes gives reliable results and second order convergence under 
certain assumptions, but is not asymptotic-preserving. Thus, if second order accuracy in 
time is desired, the L-stable Runge-Kutta method has to be applied. All three methods 
are exposed to numerical tests and compared in Section |4j 

3.2.1. Implicit Euler scheme. 

Introducing the forms 

(9,x):= / & X dx, (3.26) 

Jn 

a N (*,6,x):= /A|,vI/ 5 / 2 V||0-Vj|X^, (3.27) 

Jn 



a||(6,x) := / A[|V||e-V[|xdr, a ± (Q, X ) := / A ± V±Q ■ V±xdx, (3.28) 

Jn Jn 

allows us to write the first order, implicit Euler method in the compact notation: Find 
K +1 ,q% +1 )eV h x£ h , solution of 



, {ul + \v h ) + r{a^ul + \v h ) + n {ql + \v h )+ 1 j Vi ut +1 v h ds)={u^v h ) 
{Eap) 

a \\m K, < +1 ,w h )- eo|| {ql +l ,w h ) =0, 

(3.29) 

where the non linear term (m^ +1 ) 5 ^ 2 was replaced by a first order approximation in r : 

K +1 ) 5/2 = K + 0(r)f 2 = «) 5 / 2 + 0(r). (3.30) 

A slightly different first order AP-scheme was introduced in [TH] for the resolution of the 
same temperature balance problem. There, the (P)-problem was firstly discretized in time 
(implicit Euler), then linearized by a fixed point mapping, and finally the AP reformulation 
applied. The numerical results obtained in [16] are similar to the present ones. 
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3.2.2. Linearized Crank-Nicolson scheme. 

To construct a scheme, which is second order in time, one can come to the idea to employ 

the Crank-Nicolson scheme: Find (v% + , q% + ) e Vj x 4, solution of 

(ul +1 ,v h )+r U(u n h +1/2 ,v h ) + a ll (ql + \v h ) +j f Tx u n h +1/2 v h ds) = (u n h ,v h ) 



n+l/2 n+1/2 



(3.31) 



a \\nA u h i u h > w h) - eO|| (q'h , w h ) = . 



As one can observe, we have to deal for each fixed n, with a nonlinear equation. In the linear 



n+l/2 n+1/2 



It, 



■ WftJ 



terms, one can set u\ = \ (u^ +1 + ufy. To linearize the term a\\ n i(u™ L 
however, we shall use the standard linear extrapolation method. In other words, the non- 
linearity in this last formula, (u^ ) 5 / 2 , will be replaced by a linearized second order 
approximation in r: 

5/2/1 \ 5/2 



LI 



"+1/2x5/2 



< + g « " < _1 ) + °( r ') 



< + o « - < _1 ) 



+ 0(r 2 ). 
(3.32) 
The resulting linear system reads finally: Find (u^ +1 ,q^ +1 ) E Vh x £^, solution of 
' (< +1 ,^) + | (a ± (ul +1 ,v h ) +lf T± uf 1 v h d S > ) +ra\ { {ql + \v h ) 



(CN 



«,v h ) 



AP 



a± (v%, v h ) + 7 J r , u h Vh ds 



\ Hnl (| (3^-^r 1 ) ,ul +1 ,w h ) -ea {l (ql +1 ,w h ) 



(3.33) 



= — |o||«j (| (3< - < : ) , <, w/,) . 

Unfortunately, this method is not Asymptotic-Preserving. For small values of e one 
expects that the solution will immediately fall into the space of functions almost constant 
in the direction of the anisotropy, no matter what initial condition was imposed. In the 
case of the Crank-Nicolson scheme for large time steps compared to ej[u r ^fl 2 ^ the second 



equation in (3.33) will constrain the numerical solution to oscillate if the initial condition 



is not already in the suitable space. This requires the restrictive choice of a time step of 
the order of e/(u^) 5//2 , which yields the method inapplicable in general cases. In other 
words, the Crank-Nicolson scheme is unable to model diffusion processes for large At, due 
to the inadequate approximation of the damping processes. The Crank-Nicolson scheme 
is A-stable but not L-stable and the AP-property of a scheme is strongly related to the 
L-stability of the scheme. 

As an example of the non-convergence of the (CN ap) scheme in a general case, we show 
some numerical results corresponding to a test case defined in the Section 4.2.2 The 
initial condition is a Gaussian peak located in the center of the computational domain 
with a maximum of 10 5 if . If the time step r is too large, than u^ will immediately reach 
negative values and thus the numerical algorithm will fail in the next iteration. However, 
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if t is sufficiently small the (CNap) scheme is of second order in time. Unfortunately, the 
biggest time step that does not provoke oscillations in the numerical solution, is of the 
order of 10~ 16 s, for an initial Gaussian peak of 10 5 K. This makes the (CN^p) scheme of 
no practical use in real simulations. These results are plotted on Figure [2j 




1.85000 
160000 



120000 



80000 



50000 



(a) t = 0.1 



(b) 



= TO" 16 



Figure 2. Non convergence of the (CNap) scheme. Negative values of u\ 

are obtained after one iteration of the method, for big time steps. If the 

time step is sufficiently small, the method converges. 



3.2.3. L-stable Runge-Kutta method. 



As we are interested in an AP-scheme, which is second order accurate in time, we propose 
now a two stage Diagonally Implicit Runge-Kutta (DIRK) second order scheme, which 
does not suffer from the limitations of the Crank-Nicolson discretization. The scheme is 
developed according to the following Butcher's diagram: 



(3.34) 



A 
1 


A 
1-A A 




1-A A 



with A = 1 



i 

n/2" 
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Remark 3.1. (Butcher's diagram) The coefficients of the s-stage Runge-Kutta method are 
usually displayed in a Butcher's diagram : 



(3.35) 



C\ 


an ■ 


■ ais 


c s 


a si ■ 


&SS 




h ■ 


■ b s 



Applying this method to approximate to following problem 

du 



dt 



Lu + f{t) 



(3.36) 



reads: For given u n , being an approximation of u(t n ), the u n+1 is determined accordingly 
to : 



Ui = u 



+ T Y1 a iA Lu i + f(t + C J T ))> 



i=i 



u n+1 = u n 






j u j- 



(3.37) 
(3.38) 



If bj = a s j for j — 1, . . . , s than u 



n+l 



Us 



The scheme (3.34) is known to be L-stable, thus providing the Asymptotic Preserving 
property. The scheme writes: Find (u 1 ^ ,q% + ) 6V/,x C h , solution of 

' Kt 1 ; v h) + r\ (o ± K+ x , v h ) + 7 J r± <JV ds + a\\ (q§\v h ) 

= «,v h ) 
{ on,,, K + AK-<^),<+>,) -e H {q^,w h ) =0 



{RK 



AP 



' K£\ Vh) + rX (a±{u n 2 + h \ v h ) + 7 J r± u^v h ds + a„ (q^+\ v h ] 
= (ulv h ) + ^(u n 1 + 1 -u n h ,v h ) 



(3.39) 



■u 



a\\ni 



n+l 



(< 



+ {u n h -u 



n—l\ „,n+l 



U 



n+l 
2,h ' 



h u h 



€ +1 



),ulX 1 ,w h ) -ea\\(q%+\w h ) = 



n+l 
H2,h 



with u™^ 1 (respectively w^ 1 ) being the solution of the first (respectively second) stage of 
the Runge-Kutta method. The terms u^ + X(u^ — w^ _1 ) and u^ + (u^ — u^~ l ) are respectively 
the second order time- approximations of Uh(t + At) and Uhit + r), used to linearize the 
problem. 
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For each time step we have therefore to assemble and solve two linearized problems. This 
method is two times slower than the Crank-Nicolson scheme, with the advantage however 
of maintaining the AP-property of the scheme, advantage which is crucial for < e <C 1. 

4. Numerical results 

In this section we compare the proposed implicit Euler-AP and DIRK-AP schemes with a 
standard linearized implicit Euler discretization of the initial singular perturbation problem 



(2.3), given by 



(P) hT (u n h +1 ,v h ) + r(a ± (u n h +1 ,v h ) + ^a| N «X + \0 +7 J K +1 v h ds\ = (u n h ,v h ). 

(4.40) 

4.1. Discretization. Let us present the space discretization in a 2D case. We consider 
a square computational domain Q = [0,1] x [0,1]. All simulations are performed on 
structured meshes. Let us introduce the Cartesian, homogeneous grid 

Xi = i/N x , < i < N x , Vj = j/N y , 0<j<N y , (4.41) 

where N x and N y are positive even constants, corresponding to the number of discretization 
intervals in the x- resp. ^-direction. The corresponding mesh-sizes are denoted by h x > 
resp. h y > 0. Choosing a Q2 finite element method (Q2-FEM), based on the following 
quadratic base functions 

" ('-''-ffi-*'- x G [x t . 2 ,xl f <»-*-ffi-*-i> y G [y^yj], 

8*i = < ( " +2 "lC" +1 - x) xe[x i} x i+2 ], , e yj = \ iv^-v^-v) ye[ yj ,y j+2 ], 



y 

else I else 



(4.42) 



for even i,j and 

(x i+1 -x)(x-Xj-i) 




X G [Xi—ijXi+il, 

else ' % 



for odd i,j, we define the space 

W h :={v h = J2 v ij e *i ( x ) 9 Vi (V)} ■ 

The spaces Vh and Ch are then defined by 

Vh = Wh , C h = {qh e V h , such that q h \r in = 0}. 



[yj+i-y){y-yj-i) 


v e [yj- 


-1,1/i+i]. 





else 


(4.43) 
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The matrix elements are computed using the 2D Gauss quadrature formula, with 3 points 
in the x and y direction: 

-1 rl 1 

(AAA) 



/ / f(£,v) d £dr]= V UiUjffaTij) 
J-i J-i iJ= _ 1 



where £ — Vo — 0, £±i = r]±i = ±y|, w = 8/9 and u±i = 5/9, which is exact for 
polynomials of degree 5. Linear systems obtained for all methods in these numerical 
experiments are solved using a LU decomposition, implemented by the MUMPS library. 

4.2. Numerical tests. 

4.2.1. Known analytical solution. 

First, let us construct a numerical test case with a known solution. Finding an analytical 
solution for an arbitrary 6-field presents a considerable difficulty In the previous paper 
[5J, we presented a way to find such a solution. Let us recall briefly how to do this. The 
starting point is a limit solution 



u 







(cos (ny + a(y 2 - y) cos(7rx)) + 4) T m e t , (4.45) 



where a is a numerical constant aimed to control the variations of b. For a = 0, the limit 
solution represents a solution for the constant b case. The parameter T m is the scaling of 
u . 

Since u° is a limit solution, it is constant along the b field lines. Therefore we can 
determine the b field using the following implication 

„ n , du° , du° 

V||u° = =► h *-fa+ h y-Q- ={) i ( 4 - 46 ) 

which yields for example 

b = — B = ( a ^ V ~~ ^ cos ( 7rar ) + 7X \ (4.47) 

\B\ ' \ na(y 2 — y)sm(nx) J 

Note that the field B, constructed in this way, satisfies div5 = 0, which is an important 
property in the framework of plasma simulations. Furthermore, we have B 7^ in the 
computational domain. Now, we choose u e to be a function that converges, as e — > 0, to 
the limit solution u°, for example 

p = (cos (ny + a(y 2 - y) cos(tcx)) + 4) T m e~ l (4.48) 

q = p~ 3 / 2 sin(37T2;)/37r (4.49) 

u = p + eq. (4.50) 
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In our simulations we set a = 1 so that the direction of the anisotropy is variable in the 
computational domain. Note that we have 



— u 5/2 (t,0)V||u(t,0) = u(t,0) 
-m 5/2 (£,1)V||m(£,1) =u(t,l). 



(4.51) 
(4.52) 



The problem is supplied with a force term computed accordingly. 

As an initial condition we take u(t = 0), with u defined by (4.50). In this setting we 



expect both Asymptotic-Preserving methods (Eap) and (RKap) to converge in the optimal 
rate, independently on e and b. 

First we test the space convergence of the methods. To do this we choose a small time 
step such that the time discretization error is much smaller than the space discretization 
error. We then vary the mesh size and perform simulations for 100 time steps. The 
results are summarized in Table [T] and Figure [3j All three methods give as expected the 
third order space convergence in the L 2 -norm for large values of e. Moreover, due to the 
extremely small time step, the numerical precision is the same, even if one uses second or 
first order methods. For small values of e only the Asymptotic Preserving schemes give 
good numerical solutions. 



10 

1 


- ""\~ 


(P) -+- 

(E-MM) — •— 
(RK-MM) — •— " 


10 
1 




(P) — 
(E-MM) — •— 
(RK-MM) — •— " 


0.1 




0.1 




0.01 




N ^ 


0.01 




\ 


0.001 ' 






0.001 




H 


0.0001 






0.0001 




\ 


1e-05 






1s- 05 




\ 


1e-06 






1e-06 




> s 


1e-07 






1e-07 






1e-1E 18-10 


1e-05 1 


16-15 


1e-10 1e-05 1 




(a) h = 0.1 






(b)h 


= 0.00625 



Figure 3. Relative L 2 -errors between the exact solution u e and the 

computed solution for the standard scheme (P), Euler AP method (Eap) 

and DIRK AP scheme (RKap) as a function of e and for h = 0.1 resp. 

h = 0.00625. The time step is r = 10" 6 . 

Finally we test the time convergence of the methods. To do this we choose a small 
mesh size such that the space discretization error is smaller than the time discretization 
error. We then vary the time step and perform simulations on a fixed grid. The results are 
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L 2 -error e — 1 



0.1 



1.60 x 10" 



0.05 



2.02 x 10" 4 



0.025 2.55 x 10- 



0.0125 3.2 x 10~ 6 
0.00625 4.0 x 10~ 7 



Eap 



1.60 x 10" 



2.02 x 10" 4 



2.55 x 10" 5 



3.2 x 10~ 6 
4.0 x 10~ 7 



RKai 



1.60 x 10~ 3 



2.02 x 10" 4 



2.55 x 10" 5 



3.2 x 10~ 6 
4.0 x 10~ 7 



h 



L 2 -error e = 10" 



10 



0.1 



7.3 x 10 



-i 



0.05 



7.3 x 10- 



0.025 



7.3 x 10" 



0.0125 4.9 x 10 



-i 



0.00625 1.04 x 10" 



Eap 



1.47 x 10" 



2.04 x 10" 4 



2.65 x 10" 



3.3 x 10" 6 



RKai 



1.47 x 10" 



2.04 x 10" 4 



2.65 x 10" 



3.3 x 10" 6 



4.2 x 10 



-7 



4.2 x 10~ 7 

Table 1 . The absolute error of u in the L 2 -norm for different mesh sizes 

and e = 1 or e = 10 -10 , using the singular perturbation scheme (P) and the 

two proposed AP-schemes for a time step of r = 10 _6 s and at instant 

t = 10" 4 , with T m = 1. 



summarized in Table [2] and Figures 0] and [5j Note that the (RKap) scheme is of second 
order in time as long as the error due to the time discretization dominates the error induced 
by the space discretization. The standard (P)-scheme works well and is of first order, as 
long as e is close to one. The (Eap) scheme is of first order for all values of the anisotropic 
parameter. Also note that while the (RKap) scheme demands twice more computational 
time than the (Eap) scheme, it gives much better precision. In order to achieve a relative 
error of the order of 10~ 4 for e — 1 it suffices to take a time step of r = 0.05 in the RK- 
scheme. A comparable accuracy with (Eap) is obtained for a time step 16 times smaller. 
In the case of e = 10~ 10 the ratio is 32. 

To conclude, one can remark that the asymptotic-preserving schemes, (Eap) and (RKap), 
are uniformly accurate with respect to the perturbation parameter e. This essential fea- 
ture can be very useful in situations where the anisotropy is variable in space, i.e. the 
parameter e(x) is x-dependent. No mesh-adaptation is any more needed in these cases, a 
simple Cartesian grid enables accurate results, with no regard to the e-values. 



4.2.2. Initial Gaussian peak. 



AP-SCHEME FOR THE ANISOTROPIC HEAT EQUATION 



25 




1s-15 1e-10 



Figure 4. Relative L 2 -errors between the exact solution u e and the 

computed solution with the standard scheme (P), the Euler-AP method 

(Eap) and the DIRK-AP scheme (RKap) as a function of e and for 

t = 0.00625. The spacial grid is 200 x 200. 




(P) - 

(E-MM) — 
-. . . (RK-MM) — 




(b) e = 10" 



FlGURE 5. Relative L 2 -errors between the exact solution u e and the 

computed solution with the standard scheme (P), the Euler-AP method 

(Eap) and the DIRK-AP scheme (RKap) as a function of r, and for e = 1 

resp. e = 10~ 10 and a mesh with 200 x 200 points. Note that for e = 1 the 

P-scheme and the Eap scheme give the same precision. 



The second investigated test is the evolution of the following initial Gaussian peak, 
located in the middle of the computational domain: 



u(t = 0) 



l in 

Y 



l + e 



-50(x-0.5) 2 -50(i/-0.5) 2 



(4.53) 
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T 


L 2 -error e — 1 


P 


Eap 


RK AP 


0.1 


1.57 x 10" 2 


1.57 x 10" 2 


2.52 x 10" 3 


0.05 


8.28 x 10" 3 


8.28 x 10~ 3 


1.93 x 10" 4 


0.025 


4.25 x 10" 3 


4.25 x 10" 3 


2.62 x 10" 5 


0.0125 


2.37 x 10" 3 


2.37 x 10~ 3 


6.54 x 10~ 6 


0.00625 


1.08 x 10" 3 


1.08 x 10~ 3 


1.50 x 10" 6 


0.003125 


5.44 x 10" 4 


5.44 x 10" 4 


4.08 x 10~ 7 


0.0015625 


2.76 x 10" 4 


2.76 x 10" 4 


2.07 x 10" 7 


r 


L 2 -error e = 10~ 10 


P 


Eap 


RK AP 


0.1 


6.14 x 10" 1 


1.57 x 10" 2 


2.90 x 10" 4 


0.05 


6.30 x 10" 1 


8.22 x 10" 3 


7.21 x 10" 5 


0.025 


6.92 x lO" 1 


4.22 x 10~ 3 


1.80 x 10" 5 


0.0125 


7.08 x 10- 1 


2.36 x 10~ 3 


4.91 x 10" 6 


0.00625 


7.26 x lO" 1 


1.08 x 10~ 3 


1.15 x 10" 6 


0.003125 


7.42 x 10" 1 


5.40 x 10" 4 


3.43 x 10~ 7 


0.0015625 


6.42 x 10" 1 


2.74 x 10" 4 


2.05 x 10~ 7 



Table 2. The absolute error of u in the L 2 -norm for different time step 

using the singular perturbation scheme (P) and two proposed AP-schemes 

for mesh size 200 x 200 at time t = 0.1s with T m = 1. 



where T m = 10 5 K is the maximal temperature in the domain and the anisotropy direction 
is given as in the previous tests. We perform numerical experiments with the choice 
of £ — 1. We choose the time step r = 0.01 and perform numerical simulations on a 
fixed 50 x 50 grid with the final time set to 15s. The time step is big compared to 
the time scale induced by the initial condition. Indeed, after the first iteration of the 
algorithm the numerical solution immediately falls into the space of functions which are 
almost constant in the direction of the anisotropy (see Figure [7]). The evolution of the 
numerical solution consists of two phases. The first one, in which the parallel components 
of the diffusion operator dominate, is characterized by the exponential decay of ||w/i||i 2 (n), 
min(ti/j) and max(w/ l ) (see Figure pi). When Uh reaches some critical value, the parallel part 
of the diffusion operator becomes smaller than the perpendicular one. The direction of the 
strong diffusion is now inverted and the numerical solution aligns itself rather with the 
perpendicular direction. The minimum, maximum as well as the L 2 -norm of Uh continue 
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to approach zero, but the decay is no longer exponential. The L 2 -norm and the maximal 
value remain close to each other and almost constant in time. The minimal value of Uh, as 
well as the boundary-values decrease much faster. 




Figure 6. min(uh), max(w/ l ) and ||m/i||l 2 (£3) as a function of time for the 
Gaussian peak experiment, for T m = 10 5 and e — 1. Time step is r = 0.01s 

and a mesh size of 50 x 50. 



5. Conclusion 

The here presented Asymptotic-Preserving scheme proves to be an efficient, general and 
easy to implement numerical method for solving nonlinear, strongly anisotropic parabolic 
problems. This kind of problems occur in several important applications, as for example 
magnetically confined fusion plasmas. The method is based on a reformulation of the prob- 
lem, initially introduced by the authors in an elliptic framework, and a careful linearization 
as well as time-discretization of the resulting equation, which does not destroy the AP- 
properties of the space-discretization. Numerical experiments show clearly the advantages 
of such an AP-scheme. 
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(a) t = 



(b) t = 0.01 




(c) t = 4.5 



(d) i = 4.75 




(e) i = 5 



(f) t = 6 



Figure 7. Numerical solution at different time steps for the Gaussian 
peak experiment, for T m = 10 5 and e — 1. Time step is r = 0.01s and a 

mesh size of 50 x 50. 



ap-scheme for the anisotropic heat equation 29 

References 

[1] D. Aronson. The porous medium equation. A. Fasano, M. Primicerio (Eds.), Nonlinear Diffusion 

Problems, Lecture Notes in Mathematics, 1224:1-46, 1986. 
[2] S. F. Ashby, W. J. Bosl, R. D. Falgout, S. G. Smith, A. F. Tompson, and T. J. Williams. A Nu- 
merical Simulation of Groundwater Flow and Contaminant Transport on the CRAY T3D and C90 
Supercomputers. International Journal of High Performance Computing Applications, 13(l):80-93, 
1999. 
[3] P. Basscr and D. Jones. Diffusion-tensor mri: theory, experimental design and data analysis-a tech- 
nical review. NMR in Biomedicine, 15(7-8):456-467, 2002. 
[4] C. Bcaulicu. The basis of anisotropic water diffusion in the nervous system-a technical review. NMR 

in Biomedicine, 15(7-8) :435-455, 2002. 
[5] B. Berkowitz. Characterizing flow and transport in fractured geological media: A review. Advances 

in Water Resources, 25(8-12):861-884, 2002. 
[6] P. Degond, F. Deluzet, A. Lozinski, J. Narski, and C. Negulcscu. Duality-based asymptotic-preserving 

method for highly anisotropic diffusion equations. arXiv:1008.3405vl, 2010. 
[7] P. Degond, A. Lozinski, J. Narski, and C. Negulescu. An asymptotic-preserving method for highly 

anisotropic elliptic equations based on a micro-macro decomposition. arXiv:1102.0904vl, 2012. 
[8] Y. Dubinskii. Some integral inequalities and the solvability of degenerate quasi-linear elliptic systems 

of differential equations. Matematicheskii Sbornik, 106(3) :458-480, 1964. 
[9] Y. Dubinskii. Weak convergence for nonlinear elliptic and parabolic equations. Matematicheskii 
Sbornik, 109(4) :609-642, 1965. 
[10] S. Giinter, K. Lackner, and C. Tichmann. Finite element and higher order difference formulations 
for modelling heat transport in magnetised plasmas. Journal of Computational Physics, 226(2) :2306- 
2316, 2007. 
[11] H. Jian and B. Song. Solutions of the anisotropic porous medium equation in R™ under an ^-initial 

value. Nonlinear analysis, 64(9):2098-2111, 2006. 
[12] S. Jin. Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. SIAM J. 

Sci. Comput., 21(2):441-454, 1999. 
[13] J. Lions. Quelques methodes de resolution des problemes aux limites non lineaires. Gauthicr-Villars, 

1969. 
[14] J.-L. Lions and E. Magenes. Non-homogeneous boundary value problems and applications. Vol. I. 
Springer- Verlag, New York, 1972. Translated from the French by P. Kenneth, Die Grundlehren der 
mathematischen Wisscnschaftcn, Band 181. 
[15] H. Lutjens and J. Luciani. The xtor code for nonlinear 3d simulations of mhd instabilities in tokamak 

plasmas. Journal of Computational Physics, 227(14) :6944-6966, 2008. 
[16] A. Mentrelli and C. Negulcscu. Asymptotic preserving scheme for highly anisotropic, nonlinear diffu- 
sion equations, application: Sol plasmas, submitted to JCP, 2012. 
[17] W. Park, E. Belova, G. Fu, X. Tang, H. Strauss, and L. Sugiyama. Plasma simulation studies using 

multilevel physics models. Physics of Plasmas, 6:1796, 1999. 
[18] P. Perona and J. Malik. Scale-space and edge detection using anisotropic diffusion. Pattern Analysis 

and Machine Intelligence, IEEE Transactions on, 12(7):629-639, 1990. 
[19] M. Pierre. Personal communication, 2011. 

[20] J. Simon. Compact sets in the space LP(0,T;B). Ann. Mat. Pura Appl. (4), 146:65-96, 1987. 
[21] P. Tamain. Etude des flux de matiere dans le plasma de bord des tokamaks. PhD thesis, Marseille 1: 
2007., 2007. 



30 A. LOZINSKI, J. NARSKI, C. NEGULESCU 

[22] J. Vazquez. The porous medium equation: mathematical theory. Oxford University Press, USA, 2007. 
[23] J. Weickcrt. Anisotropic diffusion in image processing. European Consortium for Mathematics in 

Industry. B. G. Teubner, Stuttgart, 1998. 
[24] J. Wesson. Tokamaks. Oxford University Press, New York, NY, 1987. 

Universite de Toulouse, UPS, INSA, UT1, UTM, Institut de Mathematiques de Toulouse, 
118 route de Narbonne, F-31062 Toulouse, France 

E-mail address: alexei . lozinskiOmath.univ-toulouse . f r , jacek.narskiOmath.univ-toulouse . f r , 
claudia . negulescuOmath . univ-toulouse . f r 



